Larval respirometry (photosynthesis & respiration) analysis and plotting
AS Huffmyer
2021
This script analyzes and plots data for Symbiontic Integration 2021 respirometry data for Pacuta larvae.
Setup
Set up workspace, set options, and load required packages.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse')
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car')
if ("lme4" %in% rownames(installed.packages()) == 'FALSE') install.packages('lme4')
if ("lmerTest" %in% rownames(installed.packages()) == 'FALSE') install.packages('lmerTest')
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales')
if ("cowplot" %in% rownames(installed.packages()) == 'FALSE') install.packages('cowplot')
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2')
if ("effects" %in% rownames(installed.packages()) == 'FALSE') install.packages('effects')
if ("emmeans" %in% rownames(installed.packages()) == 'FALSE') install.packages('emmeans')
if ("multcomp" %in% rownames(installed.packages()) == 'FALSE') install.packages('multcomp')
#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('lme4')
library('lmerTest')
library('scales')
library('cowplot')
library('effects')
library('emmeans')
library('multcomp')Data manipulation
Load data from LoLinR.
PRdata<-read.csv("Pacu2021/Output/Respiration/PR/oxygen_P_R_calc.csv") #load dataFormat data columns.
#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]
#format columns
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Plate<-as.factor(PRdata$Plate)
PRdata$Run<-as.factor(PRdata$Run)Calculate a P:R ratio using gross photosynthesis and respiration calculated as GP:R.
PRdata$ratio<-abs(PRdata$GP.nmol.org.min)/abs(PRdata$R.nmol.org.min) #calculate ratio with absolute valuesView outliers. Remove outliers for P:R data.
boxplot(PRdata$R.nmol.org.min)
boxplot(PRdata$P.nmol.org.min)
boxplot(PRdata$GP.nmol.org.min)
boxplot(PRdata$ratio)
PRdata<-PRdata%>%filter(ratio < 4) #filter out outlier for PR ratio
boxplot(PRdata$ratio)
Plotting
Respiration
Plot data as a box plot.
r_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = R.nmol.org.min, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.2, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot 
PRdata$Date<-as.factor(PRdata$Date)
r_plot1<-PRdata %>%
group_by(Temperature, as.factor(Date))%>%
ggplot(., aes(x = as.factor(Temperature), y = R.nmol.org.min, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.2, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot1
Plot data as a mean dot plot.
r_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(R.nmol.org.min), sd=sd(R.nmol.org.min), N=length(R.nmol.org.min), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.1, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot_b
Photosynthesis (Net)
For photosynthesis filter out runs 7 and 8, these runs were outliers that didn’t work due to equipment issues.
Plot data as a box plot.
p_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = P.nmol.org.min, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); p_plot
p_plot1<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = P.nmol.org.min, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); p_plot1
Plot data as a mean dot plot.
p_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(P.nmol.org.min), sd=sd(P.nmol.org.min), N=length(P.nmol.org.min), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.005, 0.15), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); p_plot_b
Photosynthesis (Gross)
Plot data as a box plot.
gp_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = GP.nmol.org.min, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); gp_plot
gp_plot1<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = GP.nmol.org.min, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); gp_plot1
Plot data as a mean dot plot.
gp_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(GP.nmol.org.min), sd=sd(GP.nmol.org.min), N=length(GP.nmol.org.min), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.005, 0.15), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); gp_plot_b
P:R Ratio
Plot data as a box plot.
pr_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = ratio, group=Temperature)) +
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("P:R")))) +
scale_y_continuous(limits=c(0, 4), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); pr_plot
pr_plot1<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = ratio, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("P:R")))) +
scale_y_continuous(limits=c(0, 4), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); pr_plot1
Plot data as a mean dot plot.
pr_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(ratio), sd=sd(ratio), N=length(ratio), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("P:R")))) +
scale_y_continuous(limits=c(0, 2), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); pr_plot_b
Generate final figure
Box plot assembled figure.
full_fig<-plot_grid(r_plot, p_plot, gp_plot, pr_plot, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")
ggsave(filename="Pacu2021/Figures/Respiration/PR/respirometry_fig_boxplot.png", plot=full_fig, dpi=500, width=18, height=6, units="in")Mean plot assembled figure.
full_fig_b<-plot_grid(r_plot_b, p_plot_b, gp_plot_b, pr_plot_b, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")
ggsave(filename="Pacu2021/Figures/Respiration/PR/respirometry_fig_mean.png", plot=full_fig_b, dpi=500, width=18, height=6, units="in")Analysis
Respiration
Change data columns to factors.
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Date<-as.factor(PRdata$Date)
PRdata$Run<-as.factor(PRdata$Run)View respiration data.
hist(PRdata$R.nmol.org.min)
Build model and examine for respiration across temperatures.
Rmodel1<-aov(R.nmol.org.min~Temperature, data=PRdata)
summary(Rmodel1)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.00617 0.003087 2.537 0.082 .
## Residuals 180 0.21903 0.001217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check for effects of day.
Rmodel1b<-aov(R.nmol.org.min~Temperature*Date, data=PRdata)
summary(Rmodel1b)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.00617 0.003087 2.562 0.080 .
## Date 2 0.00454 0.002269 1.884 0.155
## Temperature:Date 1 0.00123 0.001228 1.019 0.314
## Residuals 177 0.21326 0.001205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check assumptions of model for residual normality and variance.
qqPlot(residuals(Rmodel1))
## [1] 131 123leveneTest(residuals(Rmodel1)~Temperature, data=PRdata)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 2.0868 0.1271
## 180Conduct post hoc test.
emm<-emmeans(Rmodel1, ~Temperature)
pairs(emm)## contrast estimate SE df t.ratio p.value
## 27 - 31 0.01214 0.00590 180 2.056 0.1022
## 27 - 34 -0.00127 0.00688 180 -0.184 0.9815
## 31 - 34 -0.01340 0.00743 180 -1.805 0.1708
##
## P value adjustment: tukey method for comparing a family of 3 estimatesNet Photosynthesis
Build model and examine for Photosynthesis across temperatures.
Pmodel1<-aov(P.nmol.org.min~Temperature, data=PRdata)
summary(Pmodel1)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.0157 0.007851 7.986 0.000475 ***
## Residuals 180 0.1769 0.000983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check for effect of day.
Pmodel1b<-aov(P.nmol.org.min~Temperature*Date, data=PRdata)
summary(Pmodel1b)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.01570 0.007851 8.328 0.000349 ***
## Date 2 0.01008 0.005040 5.347 0.005565 **
## Temperature:Date 1 0.00001 0.000009 0.009 0.922918
## Residuals 177 0.16685 0.000943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check assumptions of model for residual normality and variance.
qqPlot(residuals(Pmodel1))
## [1] 88 58leveneTest(residuals(Pmodel1)~Temperature, data=PRdata)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 7.2644 0.0009246 ***
## 180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conduct post hoc test.
emm<-emmeans(Pmodel1, ~Temperature)
pairs(emm)## contrast estimate SE df t.ratio p.value
## 27 - 31 0.00579 0.00531 180 1.090 0.5215
## 27 - 34 0.02466 0.00618 180 3.989 0.0003
## 31 - 34 0.01888 0.00667 180 2.828 0.0144
##
## P value adjustment: tukey method for comparing a family of 3 estimatesGross Photosynthesis
Build model and examine for gross photosynthesis across temperatures.
GPmodel1<-aov(GP.nmol.org.min~Temperature, data=PRdata)
summary(GPmodel1)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.0247 0.012359 4.241 0.0159 *
## Residuals 180 0.5246 0.002914
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check for effect of day.
GPmodel1b<-aov(GP.nmol.org.min~Temperature*Date, data=PRdata)
summary(GPmodel1b)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.0247 0.012359 4.395 0.0137 *
## Date 2 0.0253 0.012660 4.501 0.0124 *
## Temperature:Date 1 0.0014 0.001445 0.514 0.4745
## Residuals 177 0.4978 0.002812
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check assumptions of model for residual normality and variance.
qqPlot(residuals(GPmodel1))
## [1] 41 131leveneTest(residuals(GPmodel1)~Temperature, data=PRdata)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.2835 0.2796
## 180Conduct post hoc test.
emm<-emmeans(GPmodel1, ~Temperature)
pairs(emm)## contrast estimate SE df t.ratio p.value
## 27 - 31 -0.00635 0.00914 180 -0.695 0.7666
## 27 - 34 0.02593 0.01065 180 2.436 0.0417
## 31 - 34 0.03228 0.01149 180 2.809 0.0152
##
## P value adjustment: tukey method for comparing a family of 3 estimatesP:R
Build model and examine for PR ratio across temperatures.
PRmodel1<-aov(ratio~Temperature, data=PRdata)
summary(PRmodel1)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 2.89 1.4471 4.087 0.0184 *
## Residuals 180 63.73 0.3541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check for effect of day.
PRmodel1b<-aov(ratio~Temperature*Date, data=PRdata)
summary(PRmodel1b)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 2.89 1.4471 4.339 0.01447 *
## Date 2 4.62 2.3078 6.920 0.00128 **
## Temperature:Date 1 0.08 0.0837 0.251 0.61701
## Residuals 177 59.03 0.3335
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check assumptions of model for residual normality and variance.
qqPlot(residuals(PRmodel1))
## [1] 159 94leveneTest(residuals(PRmodel1)~Temperature, data=PRdata)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 4.1833 0.01676 *
## 180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conduct post hoc test.
emm<-emmeans(PRmodel1, ~Temperature)
pairs(emm)## contrast estimate SE df t.ratio p.value
## 27 - 31 0.189 0.101 180 1.873 0.1495
## 27 - 34 0.312 0.117 180 2.660 0.0230
## 31 - 34 0.124 0.127 180 0.975 0.5936
##
## P value adjustment: tukey method for comparing a family of 3 estimatesGenerate summary
Generate summary of all respiration data.
summary<-PRdata%>%
group_by(Temperature)%>%
filter(R.nmol.org.min<0)%>%
filter(P.nmol.org.min>0)%>%
filter(GP.nmol.org.min>0)%>%
summarise(N=length(R.nmol.org.min),
Mean_R=mean(R.nmol.org.min),
SE_R=sd(R.nmol.org.min)/sqrt(length(R.nmol.org.min)),
Mean_P=mean(P.nmol.org.min),
SE_P=sd(P.nmol.org.min)/sqrt(length(P.nmol.org.min)),
Mean_GP=mean(GP.nmol.org.min),
SE_GP=sd(GP.nmol.org.min)/sqrt(length(GP.nmol.org.min)),
Mean_PR=mean(ratio),
SE_PR=sd(ratio)/sqrt(length(ratio)))
summary## # A tibble: 3 × 10
## Temperature N Mean_R SE_R Mean_P SE_P Mean_GP SE_GP Mean_PR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 27 81 -0.0728 0.00345 0.0577 0.00358 0.130 0.00531 1.91
## 2 31 54 -0.0828 0.00374 0.0485 0.00367 0.131 0.00580 1.66
## 3 34 32 -0.0738 0.00690 0.0308 0.00295 0.105 0.00799 1.55
## # … with 1 more variable: SE_PR <dbl>summary%>%
write_csv(., "Pacu2021/Output/Respiration/PR/mean_respiration.csv")